Explanation of machine learning models predictions


Anton Kulesh, Data Scientist
Datathon
Minsk, July 2019


In this notebook we play with different tools for interpretatitng of machine learning models. We solve binary classification problem, try to predict probability of customer churn and explain predictions of the black-box model. As a black-box model we use gradient boosting on decision trees (XGBoost).

Local explanation

- LIME
- ELI5
- SHAP
- InterpretML

Global feature importance

model-specific (for tree-based ensembles)

- Gain
- Splits count
- Coverage   

model-agnostic

- Permutation
- mean(|shap_values|)


Ok, let's go!

In [1]:
import os
import warnings
warnings.filterwarnings("ignore")

# LIME
import lime
from lime.lime_tabular import LimeTabularExplainer
from lime import submodular_pick
#
# ELI5
import eli5
from eli5.sklearn import PermutationImportance
#
# SHAP
import shap
#
# InterpretML
import interpret
from interpret.data import ClassHistogram
from interpret.perf import ROC
from interpret.blackbox import ShapKernel, LimeTabular, MorrisSensitivity, PartialDependence
from interpret.glassbox import ExplainableBoostingClassifier, LogisticRegression
#
# Data manipulation and modeling
import numpy as np
import pandas as pd
from sklearn.dummy import DummyClassifier
from sklearn.preprocessing import OrdinalEncoder
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score, classification_report
import xgboost as xgb
from xgboost import XGBClassifier
#
# Plotting
import matplotlib.pyplot as plt
%matplotlib inline

# For plotting nice shap /graphs
shap.initjs()
/home/a_kulesh/anaconda3/envs/myenv/lib/python3.7/site-packages/sklearn/externals/six.py:31: DeprecationWarning: The module is deprecated in version 0.21 and will be removed in version 0.23 since we've dropped support for Python 2.7. Please rely on the official version of six (https://pypi.org/project/six/).
  "(https://pypi.org/project/six/).", DeprecationWarning)

1. Data

Telco Customer Churn

Data description

Each row represents a customer, each column contains customer's attributes. The data set includes information about:

  1. Customers who left within the last month – the column is called Churn (target variable)
  2. Services that each customer has signed up for – phone (PhoneService), multiple lines (MultipleLines), internet (InternetService), online security (OnlineSecurity), online backup (OnlineBackup), device protection (DeviceProtection), tech support (TechSupport), streaming TV (StreamingTV) and movies (StreamingMovies)
  3. Customer account information – how long they’ve been a customer (tenure), contract (Contract), payment method (PaymentMethod), paperless billing (PaperlessBilling), monthly charges (MonthlyCharges), and total charges (TotalCharges)
  4. Demographic info about customers – gender, age range (SeniorCitizen), and if they have partners (Partner) and dependents (Dependents)

https://www.kaggle.com/blastchar/telco-customer-churn

Data preparation

After loading the data we'll:

  • remove useless columns
  • fill missing values
  • encode binary variables
  • encode categorical variables
In [2]:
path_to_data = "./data/telco-customer-churn.zip"
data = pd.read_csv(path_to_data, compression="zip")

del data["customerID"]
data['TotalCharges'] = data['TotalCharges'].replace(" ", 0).astype('float32')
data['gender'] = data['gender'].apply(lambda x: 1 if x == "Female" else 0)

bool_columns = ['Partner', 'Dependents', 'PhoneService',
                'OnlineSecurity', 'OnlineBackup', 'DeviceProtection',
                'TechSupport', 'StreamingTV', 'StreamingMovies',
                'PaperlessBilling', 'Churn'
               ]
for col in bool_columns:
    data[col] = data[col].apply(lambda x: 1 if x == "Yes" else 0)

columns = data.columns.to_list()[:-1]
categorical_names = ['MultipleLines', 'InternetService', 'Contract', 'PaymentMethod']
categorical_columns = [columns.index(col) for col in categorical_names]
encoder = OrdinalEncoder()
data[categorical_names] = encoder.fit_transform(data[categorical_names])
categorical_names_dict = {col: list(values) for col, values in zip(categorical_columns, encoder.categories_)}
In [3]:
def categorical_decode(data, idx):
    return dict(zip(
            categorical_names,
            encoder.inverse_transform(
            data[idx, categorical_columns].reshape(1, -1))[0]
            )
        )
In [4]:
data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 7043 entries, 0 to 7042
Data columns (total 20 columns):
gender              7043 non-null int64
SeniorCitizen       7043 non-null int64
Partner             7043 non-null int64
Dependents          7043 non-null int64
tenure              7043 non-null int64
PhoneService        7043 non-null int64
MultipleLines       7043 non-null float64
InternetService     7043 non-null float64
OnlineSecurity      7043 non-null int64
OnlineBackup        7043 non-null int64
DeviceProtection    7043 non-null int64
TechSupport         7043 non-null int64
StreamingTV         7043 non-null int64
StreamingMovies     7043 non-null int64
Contract            7043 non-null float64
PaperlessBilling    7043 non-null int64
PaymentMethod       7043 non-null float64
MonthlyCharges      7043 non-null float64
TotalCharges        7043 non-null float32
Churn               7043 non-null int64
dtypes: float32(1), float64(5), int64(14)
memory usage: 1.0 MB

Let's look at the prepared data sample

In [5]:
data.head().T
Out[5]:
0 1 2 3 4
gender 1.00 0.00 0.000000 0.00 1.000000
SeniorCitizen 0.00 0.00 0.000000 0.00 0.000000
Partner 1.00 0.00 0.000000 0.00 0.000000
Dependents 0.00 0.00 0.000000 0.00 0.000000
tenure 1.00 34.00 2.000000 45.00 2.000000
PhoneService 0.00 1.00 1.000000 0.00 1.000000
MultipleLines 1.00 0.00 0.000000 1.00 0.000000
InternetService 0.00 0.00 0.000000 0.00 1.000000
OnlineSecurity 0.00 1.00 1.000000 1.00 0.000000
OnlineBackup 1.00 0.00 1.000000 0.00 0.000000
DeviceProtection 0.00 1.00 0.000000 1.00 0.000000
TechSupport 0.00 0.00 0.000000 1.00 0.000000
StreamingTV 0.00 0.00 0.000000 0.00 0.000000
StreamingMovies 0.00 0.00 0.000000 0.00 0.000000
Contract 0.00 1.00 0.000000 1.00 0.000000
PaperlessBilling 1.00 0.00 1.000000 0.00 1.000000
PaymentMethod 2.00 3.00 3.000000 0.00 2.000000
MonthlyCharges 29.85 56.95 53.850000 42.30 70.700000
TotalCharges 29.85 1889.50 108.150002 1840.75 151.649994
Churn 0.00 0.00 1.000000 0.00 1.000000

Let's look at the target variable distribution

In [6]:
target = "Churn"
data[target].plot(kind="hist");

We can observe some class imbalance. But it's not big problem.

In [7]:
data[target].value_counts(1)
Out[7]:
0    0.73463
1    0.26537
Name: Churn, dtype: float64

We omit feature engineering phase, as this is redundant in context of our goals.

Make dataset

Let's simply split out data to train (70% for model training) and test (30% for model validation) sets.

In [8]:
X_data = data.copy()
y_data = X_data.pop(target)
In [9]:
X_train, X_test, y_train, y_test = train_test_split(X_data.values, y_data.values, test_size=0.3, random_state=42)
print("Train size: {}".format(X_train.shape))
print("Test size: {}".format(X_test.shape))
Train size: (4930, 19)
Test size: (2113, 19)

2. Training

For modeling we use simple and powerful gradient boosting classifier. This is our black box model that we want to explain. Let's train our model on given dataset and calculate quality metric (ROC AUC).

In [10]:
clf = XGBClassifier(max_depth=3, n_estimators=50, importance_type="gain")
In [11]:
clf.fit(X_train, y_train)
y_proba = clf.predict_proba(X_test)[:, 1]
print("ROC-AUC: {}".format(roc_auc_score(y_test, y_proba)))
ROC-AUC: 0.854261896837849

Single Tree

Let's take a closer look at the first tree from our ensemble.

In [12]:
booster = clf.get_booster()
original_feature_names = booster.feature_names
booster.feature_names = columns
print(booster.get_dump()[0])
booster.feature_names = original_feature_names
0:[Contract<0.5] yes=1,no=2,missing=1
	1:[MonthlyCharges<68.0250015] yes=3,no=4,missing=3
		3:[tenure<3.5] yes=7,no=8,missing=7
			7:leaf=-0.0190954786
			8:leaf=-0.114973262
		4:[tenure<13.5] yes=9,no=10,missing=9
			9:leaf=0.0739255026
			10:leaf=-0.0477958247
	2:[MonthlyCharges<93.6750031] yes=5,no=6,missing=5
		5:[Contract<1.5] yes=11,no=12,missing=11
			11:leaf=-0.169309467
			12:leaf=-0.194810808
		6:[Contract<1.5] yes=13,no=14,missing=13
			13:leaf=-0.0985714272
			14:leaf=-0.16579926

Model performance results

We put our predictions and true labels in one table. Then we convert predicted probabilities into binary labels (threshold=0.5)

In [13]:
results = pd.DataFrame(np.vstack((y_test, y_proba)).T, columns=["y_test", "y_proba"])
results["y_test"] = results["y_test"].astype("bool")
results["y_pred"] = results["y_proba"] >= 0.5
results["error"] = results['y_proba'] - results['y_test']
results["abs_error"] = abs(results['y_proba'] - results['y_test'])
results.sort_values('error', ascending=False, inplace=True)

Let's at the classification report

In [14]:
print(classification_report(results["y_test"], results["y_pred"]))
              precision    recall  f1-score   support

       False       0.83      0.91      0.87      1539
        True       0.68      0.50      0.58       574

    accuracy                           0.80      2113
   macro avg       0.76      0.71      0.72      2113
weighted avg       0.79      0.80      0.79      2113

Now let's compare performance of our "black box" with "dummy" model (generates predictions by respecting the training set's class distribution)

In [15]:
dummy = DummyClassifier()
dummy.fit(X_train, y_train)
y_proba_dummy = dummy.predict_proba(X_test)[:, 1]
print("ROC-AUC: {}".format(roc_auc_score(y_test, y_proba_dummy)))
print(classification_report(y_test, y_proba_dummy >= 0.5))
ROC-AUC: 0.4996292673870765
              precision    recall  f1-score   support

           0       0.73      0.74      0.73      1539
           1       0.27      0.26      0.27       574

    accuracy                           0.61      2113
   macro avg       0.50      0.50      0.50      2113
weighted avg       0.60      0.61      0.61      2113

Not bad! We can see that model shows quite promising results even without applying smart feature engineering technique and model parameters tuning. So let's try to give some explanations and crack our black box open.


Local Interpretability

Most part of this notebook is devoted to the explanation of individual samples (local interpretability of the model).

We start with the LIME algorithm (original implementation)

lime

lime package provides working with models which have different input types:

  • LimeTabularExplainer -- working with tabular data
  • RecurrentTabularExplainer -- working with time series
  • LimeTextExplainer -- working with text data inputs
  • LimeImageExplainer -- explains predictions on images

For our tast we use LimeTabularExplainer. Let's create explainer (LimeTabularExplainer class object) and get explanation for some examples (explainer.explain_instance).

In [16]:
def create_explainer(data, **kwargs):
    explainer = LimeTabularExplainer(
        data, feature_names=columns, categorical_features=categorical_columns, 
        categorical_names=categorical_names_dict, class_names=["no churn", "churn"],
        discretize_continuous=True, **kwargs)

    return explainer

Example of "no churn"

Let's look at the top 5 features and their contributions

In [17]:
i = 2
explainer = create_explainer(X_train, verbose=True, kernel_width=None)

exp = explainer.explain_instance(X_test[i], predict_fn=clf.predict_proba, num_samples=1000,
                                 num_features=10, model_regressor=None)
print("Predicted_label: [{}]".format(int(y_proba[i] >= 0.5)))
print("True_label: [{}]".format(y_test[i]))
Intercept 0.29249616285592495
Prediction_local [-0.08252804]
Right: 0.016522923
Predicted_label: [0]
True_label: [0]
In [18]:
exp_df = pd.DataFrame(exp.as_list(), columns=['feature_value', 'feature_contribution'])
exp_df
Out[18]:
feature_value feature_contribution
0 Contract=Two year -0.194552
1 InternetService=No -0.097012
2 MonthlyCharges <= 35.40 -0.068921
3 OnlineSecurity <= 0.00 0.042957
4 SeniorCitizen <= 0.00 -0.040116
5 29.00 < tenure <= 55.00 -0.035540
6 TechSupport <= 0.00 0.018930
7 400.50 < TotalCharges <= 1402.67 -0.014432
8 0.00 < gender <= 1.00 0.013662
9 PhoneService <= 1.00 0.000000

Let's sum weights and compare with Prediction_local above. We can see that it's exactly the same value.

In [19]:
exp_df.feature_contribution.sum() + exp.intercept[1]
Out[19]:
-0.08252804207653852

Detailed explanation

In [20]:
exp.show_in_notebook(show_table=True, show_all=False)

Feature contributions

In [21]:
exp.as_pyplot_figure();

Descriptive statistics of numerical features

In [22]:
num_columns = ['tenure', 'MonthlyCharges', 'TotalCharges']
X_data[num_columns].describe()
Out[22]:
tenure MonthlyCharges TotalCharges
count 7043.000000 7043.000000 7043.000000
mean 32.371149 64.761692 2279.734375
std 24.559481 30.090047 2266.794434
min 0.000000 18.250000 0.000000
25% 9.000000 35.500000 398.549988
50% 29.000000 70.350000 1394.550049
75% 55.000000 89.850000 3786.599976
max 72.000000 118.750000 8684.799805

We can also save our explanation as html file

exp.save_to_file("./data/demo_0.html)

Example of "churn"

In [23]:
i = 0
exp = explainer.explain_instance(X_test[i], predict_fn=clf.predict_proba, num_samples=1000,
                                 num_features=10, model_regressor=None)
print("Predicted_label: [{}]".format(int(y_proba[i] >= 0.5)))
print("True_label: [{}]".format(y_test[i]))
Intercept -0.020225115118940956
Prediction_local [0.4461491]
Right: 0.6259035
Predicted_label: [1]
True_label: [1]
In [24]:
exp.as_pyplot_figure();

Submodular Pick

This technique allows us to select a set of representative instances from our dataset.

In [25]:
explainer = create_explainer(X_train)
%time sp_explanations = submodular_pick.SubmodularPick(explainer, X_test, clf.predict_proba, method="full")
CPU times: user 14min 40s, sys: 6min 45s, total: 21min 26s
Wall time: 8min 47s
In [26]:
for i, explanation in zip(sp_explanations.V, sp_explanations.sp_explanations):
    print("\nExplanation #{}".format(i))
    print("\nIntercept: {}".format(explanation.intercept))
    print("\nLocal_prediction: {}".format(explanation.local_pred[0]))
    print("\nModel_prediction: {}".format(y_proba[i]))
    print("\nTrue_label: {}".format(y_test[i]))
    explanation.show_in_notebook();
Explanation #2104

Intercept: {1: -0.044806852684115855}

Local_prediction: 0.5600895333399646

Model_prediction: 0.6891471147537231

True_label: 1
Explanation #1572

Intercept: {0: 0.6386623985039668}

Local_prediction: 1.1222461887883766

Model_prediction: 0.023438628762960434

True_label: 0
Explanation #88

Intercept: {0: 0.7162888738121403}

Local_prediction: 1.04454338332582

Model_prediction: 0.05333323031663895

True_label: 0
Explanation #962

Intercept: {1: -0.013099823438362218}

Local_prediction: 0.4744682056799112

Model_prediction: 0.6131873726844788

True_label: 0
Explanation #1620

Intercept: {0: 0.6942977182067483}

Local_prediction: 1.037254565184352

Model_prediction: 0.028413619846105576

True_label: 0

We can observe that for some examples local prediction has high residuals. It looks like a weakness of LIME approach. Maybe it related with neibourhoods definition and local kernel settings. So, in general, we can try to pick optimal parameters for our explainer, but it out-of-scope this notebook.

See also

ELI5

Next tool is ELI5.

Explanation of a single prediction

In [25]:
i = 0
eli5.show_prediction(clf, X_test[i],
                     feature_names=columns,
                     show_feature_values=True,
                     show=["targets",
                           "feature_importances",
                           "description",
                           "method"])
Out[25]:

y=1 (probability 0.626, score 0.515) top features

Contribution? Feature Value
+0.702 Contract 0.000
+0.583 tenure 1.000
+0.292 TotalCharges 24.800
+0.111 MultipleLines 1.000
+0.107 OnlineSecurity 0.000
+0.105 PaymentMethod 2.000
+0.093 PaperlessBilling 1.000
+0.043 TechSupport 0.000
+0.023 InternetService 0.000
+0.019 OnlineBackup 0.000
-0.034 SeniorCitizen 0.000
-0.040 StreamingMovies 0.000
-0.360 MonthlyCharges 24.800
-1.130 <BIAS> 1.000
Features with largest coefficients.

Feature weights are calculated by following decision paths in trees
of an ensemble. Each leaf has an output score, and expected scores can also be
assigned to parent nodes. Contribution of one feature on the decision path
is how much expected score changes from parent to child. Weights of all 
features sum to the output score of the estimator.

Caveats:
1. Feature weights just show if the feature contributed positively or
   negatively to the final score, and does not show how increasing or
   decreasing the feature value will change the prediction.
2. In some cases, feature weight can be close to zero for an important feature.
   For example, in a single tree that computes XOR function, the feature at the
   top of the tree will have zero weight because expected scores for both
   branches are equal, so decision at the top feature does not change the
   expected score. For an ensemble predicting XOR functions it might not be
   a problem, but it is not reliable if most trees happen to choose the same
   feature at the top.

Explained as: decision paths

In [26]:
eli5.explain_prediction_df(clf, X_test[i], feature_names=columns)
Out[26]:
target feature weight value
0 1 Contract 0.702181 0.000000
1 1 tenure 0.582941 1.000000
2 1 TotalCharges 0.292155 24.799999
3 1 MultipleLines 0.111396 1.000000
4 1 OnlineSecurity 0.106951 0.000000
5 1 PaymentMethod 0.105294 2.000000
6 1 PaperlessBilling 0.093165 1.000000
7 1 TechSupport 0.042601 0.000000
8 1 InternetService 0.023000 0.000000
9 1 OnlineBackup 0.018951 0.000000
10 1 SeniorCitizen -0.034086 0.000000
11 1 StreamingMovies -0.040392 0.000000
12 1 MonthlyCharges -0.359828 24.800000
13 1 <BIAS> -1.129647 1.000000

Feature Importance

In [27]:
eli5.show_weights(clf, feature_names=columns)
Out[27]:
Weight Feature
0.3044 Contract
0.0941 InternetService
0.0931 tenure
0.0831 OnlineSecurity
0.0655 PaperlessBilling
0.0533 StreamingMovies
0.0523 MonthlyCharges
0.0507 SeniorCitizen
0.0457 TechSupport
0.0404 MultipleLines
0.0306 OnlineBackup
0.0305 TotalCharges
0.0303 PaymentMethod
0.0205 gender
0.0055 Partner
0 DeviceProtection
0 PhoneService
0 StreamingTV
0 Dependents

It's just a nice way to show model's feature importances. XGBoost by default uses gain, that is the average gain of the feature.

In [28]:
pd.DataFrame.from_records(
    data=zip(columns, clf.feature_importances_),
    columns=["name", "importance"]).sort_values(by="importance", ascending=False)
Out[28]:
name importance
14 Contract 0.304384
7 InternetService 0.094096
4 tenure 0.093078
8 OnlineSecurity 0.083126
15 PaperlessBilling 0.065494
13 StreamingMovies 0.053310
17 MonthlyCharges 0.052273
1 SeniorCitizen 0.050693
11 TechSupport 0.045686
6 MultipleLines 0.040430
9 OnlineBackup 0.030617
18 TotalCharges 0.030456
16 PaymentMethod 0.030286
0 gender 0.020528
2 Partner 0.005543
12 StreamingTV 0.000000
10 DeviceProtection 0.000000
5 PhoneService 0.000000
3 Dependents 0.000000

ELI5 also allows us to use Permutation Importance technique, and we will mention it in "Global Interpretability" section.

SHAP

SHAP (SHapley Additive exPlanations) is a unified approach to explain the output of any machine learning model. SHAP connects game theory with local explanations, uniting several previous methods and representing the only possible consistent and locally accurate additive feature attribution method based on expectations.

see repo for details: https://github.com/slundberg/shap

TreeSHAP

Efficient implementation of SHAP for tree models: https://arxiv.org/pdf/1802.03888.pdf.

Let's create explainer and calculate SHAP values

In [29]:
explainer = shap.TreeExplainer(clf, X_train, model_output="probability", feature_dependence="independent")
%time shap_values = explainer.shap_values(X_test[:500])
100%|===================| 499/500 [00:22<00:00]        
CPU times: user 22 s, sys: 60 ms, total: 22.1 s
Wall time: 22 s

SHAP values of the first sample

In [30]:
shap_values[0]
Out[30]:
array([-1.87082101e-05, -6.78353330e-03,  3.22998241e-04,  0.00000000e+00,
        1.35792022e-01,  0.00000000e+00,  1.19316652e-02,  1.97713521e-02,
        2.10519546e-02,  3.79110347e-03,  0.00000000e+00,  7.96801089e-03,
        0.00000000e+00, -1.08800055e-02,  1.55867296e-01,  2.30644137e-02,
        1.85767613e-02, -7.77787246e-02,  5.89715607e-02])
In [31]:
print("Model prediction: {}".format(y_proba[0]))
print("Sum of SHAP values + base values: {}".format(sum(shap_values[0]) + explainer.expected_value))
Model prediction: 0.6259034872055054
Sum of SHAP values + base values: 0.6259034935659924
In [32]:
print("Expected value: {}".format(explainer.expected_value))
print("Average prediction: {}".format(clf.predict_proba(X_train)[:, 1].mean(0)))
Expected value: 0.2642553270057042
Average prediction: 0.2642553150653839

Explanation of a single prediction

In [33]:
i = 0
shap.force_plot(base_value=explainer.expected_value, shap_values=shap_values[i,:], features=X_test[i,:], feature_names=columns)
Out[33]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Explanation of a multiple predictions

In [34]:
shap.force_plot(explainer.expected_value, shap_values[:500,:], X_test[:500,:], feature_names=columns)
Out[34]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

Summary Plot

[..] To get an overview of which features are most important for a model we can plot the SHAP values of every feature for every sample. The plot below sorts features by the sum of SHAP value magnitudes over all samples, and uses SHAP values to show the distribution of the impacts each feature has on the model output. The color represents the feature value (red high, blue low).
In [35]:
shap.summary_plot(shap_values, X_test[:500], feature_names=columns)

Summary Plot for "Churn" class

In [36]:
churn_mask = y_test[:500] == 1
no_churn_mask = y_test[:500] == 0
print("Customers in churn: {}".format(sum(churn_mask)))
print("Customers is alive: {}".format(sum(no_churn_mask)))
Customers in churn: 141
Customers is alive: 359
In [37]:
shap.summary_plot(shap_values[churn_mask], X_test[:500][churn_mask], feature_names=columns)

Summary plot for "No Churn" class

In [38]:
shap.summary_plot(shap_values[no_churn_mask], X_test[:500][no_churn_mask], feature_names=columns)

Feature Importances

In [39]:
shap.summary_plot(shap_values, X_test, plot_type="bar", feature_names=columns)

Dependence Plots

[..] SHAP dependence plots show the effect of a single feature across the whole dataset. They plot a feature's value vs. the SHAP value of that feature across many samples. SHAP dependence plots are similar to partial dependence plots, but account for the interaction effects present in the features, and are only defined in regions of the input space supported by data. The vertical dispersion of SHAP values at a single feature value is driven by interaction effects, and another feature is chosen for coloring to highlight possible interactions.

Below we create dependence plots for most important 7 features

In [40]:
feature_importance = pd.DataFrame.from_records(
    data=list(zip(columns, np.mean(np.abs(shap_values), axis=0))), 
    columns=["name", "importance"]
)
feature_importance.sort_values(by="importance", ascending=False, inplace=True)

n_top = 7
for col in feature_importance.name[:n_top]:
    shap.dependence_plot(col, shap_values, X_test[:500], feature_names=columns)

KernelSHAP

[..] Kernel SHAP is a method that uses a special weighted linear regression
to compute the importance of each feature. The computed importance values
are Shapley values from game theory and also coefficents from a local linear
regression.
In [41]:
model = lambda x: clf.predict_proba(x)[:, 1]

Use apply kmeans to reduce input dimensions, get approximation of original data and use it as a background dataset.

In [42]:
data = shap.kmeans(X_train, 100)
kernel_explainer = shap.KernelExplainer(model, data, feature_names=columns)
In [43]:
print("Expected values (KernalExplainer): {}".format(kernel_explainer.expected_value))
print("Expected values (TreeExplainer): {}".format(explainer.expected_value))
Expected values (KernalExplainer): 0.26739720949718
Expected values (TreeExplainer): 0.2642553270057042
In [44]:
%time shap_values_kernel = kernel_explainer.shap_values(X_test[0])
CPU times: user 1.55 s, sys: 51.9 ms, total: 1.6 s
Wall time: 1.56 s

KernelExplainer values

In [45]:
shap_values_kernel
Out[45]:
array([ 0.        ,  0.        ,  0.        ,  0.        ,  0.13395553,
        0.        ,  0.00566718,  0.00491195,  0.0113764 ,  0.00290951,
        0.        ,  0.00618924,  0.        , -0.01400079,  0.19192476,
        0.01914348,  0.01382495, -0.07308912,  0.05569319])
In [46]:
shap.force_plot(base_value=float(kernel_explainer.expected_value),
                shap_values=shap_values_kernel,
                feature_names=columns)
Out[46]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

TreeExplainer values

In [47]:
shap_values[0, :]
Out[47]:
array([-1.87082101e-05, -6.78353330e-03,  3.22998241e-04,  0.00000000e+00,
        1.35792022e-01,  0.00000000e+00,  1.19316652e-02,  1.97713521e-02,
        2.10519546e-02,  3.79110347e-03,  0.00000000e+00,  7.96801089e-03,
        0.00000000e+00, -1.08800055e-02,  1.55867296e-01,  2.30644137e-02,
        1.85767613e-02, -7.77787246e-02,  5.89715607e-02])
In [48]:
shap.force_plot(base_value=explainer.expected_value,
                shap_values=shap_values[0],
                feature_names=columns)
Out[48]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
In [49]:
mean_squared_error(shap_values[0, :], shap_values_kernel)
Out[49]:
9.408791331151606e-05

InterpretML

Microsoft Research has developed an algorithm called the Explainable Boosting Machine (EBM) which has both high accuracy and intelligibility. EBM uses modern machine learning techniques like bagging and boosting to breathe new life into traditional GAMs (Generalized Additive Models). This makes them as accurate as random forests and gradient boosted trees, and also enhances their intelligibility and editability.

This tool looks very promising. It has nice functionality but is still under development.

see repo: https://github.com/microsoft/interpret

Data Visualization

In [50]:
hist = ClassHistogram(feature_names=columns).explain_data(X_train, y_train, name="EDA")
interpret.show(hist)

Train "White Box" Model

Exaplainable Boosting Machine

Link to paper

In [51]:
ebc = ExplainableBoostingClassifier(feature_names=columns, random_state=42)
%time ebc.fit(X_train, y_train)
CPU times: user 341 ms, sys: 77 ms, total: 418 ms
Wall time: 11.2 s
Out[51]:
ExplainableBoostingClassifier(data_n_episodes=2000,
                              early_stopping_run_length=50,
                              early_stopping_tolerance=1e-05,
                              feature_names=['gender', 'SeniorCitizen',
                                             'Partner', 'Dependents', 'tenure',
                                             'PhoneService', 'MultipleLines',
                                             'InternetService',
                                             'OnlineSecurity', 'OnlineBackup',
                                             'DeviceProtection', 'TechSupport',
                                             'StreamingTV', 'StreamingMovies',
                                             'Contract', 'Paperless...
                                             'categorical', 'categorical',
                                             'categorical', 'categorical',
                                             'categorical', 'categorical',
                                             'continuous', 'categorical',
                                             'continuous', 'continuous',
                                             'continuous'],
                              holdout_size=0.15, holdout_split=0.15,
                              interactions=0, learning_rate=0.01,
                              max_tree_splits=2, min_cases_for_splits=2,
                              n_estimators=16, n_jobs=-2, random_state=42,
                              schema=None, scoring=None,
                              training_step_episodes=1)

Explanation of a single prediction

In [52]:
n = 10
ebc_local = ebc.explain_local(X_test[:n], y_test[:n], name="EBC")
interpret.show(ebc_local)
In [53]:
i = 0
categorical_decode(X_test, i)
Out[53]:
{'MultipleLines': 'No phone service',
 'InternetService': 'DSL',
 'Contract': 'Month-to-month',
 'PaymentMethod': 'Electronic check'}

Feature Importances

In [54]:
ebc_global = ebc.explain_global(name="EBC")
interpret.show(ebc_global)

Model Performance

In [55]:
ebc_roc_curve = interpret.perf.ROC(
    predict_fn=ebc.predict_proba).explain_perf(X_test, y_test, name='Churn detection (EBC)')
interpret.show(ebc_roc_curve)

All in one place

Add one more model

In [56]:
lr = LogisticRegression(feature_names=columns)
lr.fit(X_train, y_train)
lr_global = lr.explain_global(name="LR")
lr_local = lr.explain_local(X_test[:n], y_test[:n], name="LR")
lr_roc_curve = interpret.perf.ROC(
    predict_fn=lr.predict_proba).explain_perf(X_test, y_test, name='Churn detection (LR)')

Let's look at the entire picture

In [57]:
interpret.show(
    explanation=[hist, ebc_roc_curve, ebc_global, ebc_local, lr_roc_curve, lr_global, lr_local],
    share_tables=True
)
/home/a_kulesh/anaconda3/envs/myenv/lib/python3.7/site-packages/plotly/tools.py:465: DeprecationWarning:

plotly.tools.make_subplots is deprecated, please use plotly.subplots.make_subplots instead

LIME and SHAP support

Explanation of tabular data using LIME (LimeTabular)

In [58]:
n = 10
lime_explainer = LimeTabular(predict_fn=clf.predict_proba,
                             data=X_train,
                             feature_names=columns,
                             random_state=42
                            )
%time lime_local = lime_explainer.explain_local(X_test[:n], y_test[:n], name='LIME (LimeTabular)')
CPU times: user 4.39 s, sys: 1.63 s, total: 6.02 s
Wall time: 2.64 s

Currently interpret supports only ShapKernel (aka KernelShap)

In [59]:
background_data = shap.kmeans(X_train, 100).data
In [60]:
shap_explainer = ShapKernel(predict_fn=clf.predict_proba, data=background_data, feature_names=columns)
shap_local = shap_explainer.explain_local(X_test[:n], y_test[:n], name='SHAP (ShapKernel)')

In [61]:
interpret.show([lime_local, shap_local], share_tables=True)

Global Interpretability

Feature Importance techniques

Here we shortly look at some technique for assessing global feature importance.

First, let's define some helpful functions:

In [80]:
def feature_importances_table(feature_importances_array, feature_names=None, n_top=None):
    if feature_names is None:
        feature_names = X_data.columns
    if n_top is None:
        n_top = len(feature_importances_array)
    feature_importances = list(zip(feature_names, feature_importances_array))
    feature_importances = pd.DataFrame.from_records(
        feature_importances, columns=["feature_name", "importance"])
    feature_importances.sort_values("importance", inplace=True, ascending=True)
    feature_importances = feature_importances[-n_top:]
    return feature_importances


def show_feature_importances(feature_importances, feature_names=None, n_top=None):
    if not isinstance(feature_importances, pd.DataFrame):
        feature_importances = feature_importances_table(
            feature_importances, feature_names=feature_names, n_top=n_top)
        
    ax = feature_importances.plot.barh(
        x="feature_name", y="importance", figsize=(10, 8),
        fontsize=16, color='green', alpha=0.8)
    ax.grid()
    plt.title("Feature importances", fontsize=20)
    plt.show()
    

def calculate_feature_saturation(feature_importances, model=None,
                                 feature_names=None, n_top=None, show=True, label=None):
    if model is None:
        model = XGBClassifier(max_depth=3, n_estimators=50)
    if not isinstance(feature_importances, pd.DataFrame):
        feature_importances = feature_importances_table(
            feature_importances, feature_names=feature_names, n_top=n_top)
    
    score = {}
    for i in range(1, len(feature_importances)):
        features = feature_importances[-i:]["feature_name"]
        column_mask = X_data.columns.isin(features)
        model.fit(X_train[:, column_mask], y_train)
        y_proba = model.predict_proba(X_test[:, column_mask])[:, 1]
        score[i] = roc_auc_score(y_test, y_proba)
    return score


def show_feature_saturation(scores, labels, colors, show=True):
    plt.figure(figsize=(12, 6))
    for score, label, color in zip(scores, labels, colors):
        x = list(score.keys())
        y = list(score.values())
        plt.plot(x, y, "-*", color=color, label=label)
        plt.title("Dependence of the model score on the number of features")
        plt.xlabel("Numbel of features")
        plt.ylabel("ROC AUC")
    plt.legend()
    plt.grid()
    plt.show()     

Permutation Importances

In [63]:
perm = PermutationImportance(clf, random_state=42).fit(X_train, y_train)
eli5.show_weights(perm, feature_names=X_data.columns.to_list())
Out[63]:
Weight Feature
0.0551 ± 0.0036 tenure
0.0349 ± 0.0014 Contract
0.0348 ± 0.0069 MonthlyCharges
0.0175 ± 0.0028 InternetService
0.0076 ± 0.0039 OnlineSecurity
0.0070 ± 0.0024 TotalCharges
0.0070 ± 0.0017 PaperlessBilling
0.0062 ± 0.0034 PaymentMethod
0.0033 ± 0.0014 MultipleLines
0.0032 ± 0.0018 SeniorCitizen
0.0023 ± 0.0020 TechSupport
0.0008 ± 0.0010 StreamingMovies
0.0000 ± 0.0002 Partner
0 ± 0.0000 Dependents
0 ± 0.0000 gender
0 ± 0.0000 PhoneService
0 ± 0.0000 DeviceProtection
0 ± 0.0000 StreamingTV
-0.0004 ± 0.0024 OnlineBackup
In [64]:
feature_importances_by_permutation = np.abs(perm.feature_importances_)
show_feature_importances(feature_importances_by_permutation)

SHAP Values (Tree SHAP)

In [65]:
explainer = shap.TreeExplainer(clf, X_train, model_output="probability", feature_dependence="independent")
%time shap_values = explainer.shap_values(X_train)
100%|===================| 4927/4930 [03:36<00:00]        
CPU times: user 3min 35s, sys: 568 ms, total: 3min 36s
Wall time: 3min 35s
In [66]:
feature_importances_by_shap = np.mean(np.abs(shap_values), axis=0)
show_feature_importances(feature_importances_by_shap, n_top=None)

Morris method

See for details: https://en.wikipedia.org/wiki/Morris_method

In [67]:
sensitivity = MorrisSensitivity(predict_fn=clf.predict_proba, data=X_train, feature_names=columns)
sensitivity_global = sensitivity.explain_global(name="Global Sensitivity")
In [68]:
feature_importances_by_moris = sensitivity_global.data()["scores"]
show_feature_importances(feature_importances_by_moris)

Gain

The average training loss reduction gained when using a feature for splitting.

In [69]:
clf = XGBClassifier(max_depth=3, n_estimators=50, importance_type="gain")
clf.fit(X_train, y_train)

feature_importances_by_gain = clf.feature_importances_
show_feature_importances(feature_importances_by_gain, n_top=None)

Weight

The number of times a feature is used to split the data across all trees.

In [70]:
clf = XGBClassifier(max_depth=3, n_estimators=50, importance_type="weight")
clf.fit(X_train, y_train)

feature_importances_by_weight = clf.feature_importances_
show_feature_importances(feature_importances_by_weight, n_top=None)

Cover

The number of times a feature is used to split the data across all trees weighted by the number of training data points that go through those splits.

In [71]:
clf = XGBClassifier(max_depth=3, n_estimators=50, importance_type="cover")
clf.fit(X_train, y_train)

feature_importances_by_cover = clf.feature_importances_
show_feature_importances(feature_importances_by_cover, n_top=None)
In [72]:
scores = [calculate_feature_saturation(fi)
          for fi in [
              feature_importances_by_permutation,
              feature_importances_by_shap,
              feature_importances_by_moris,
              feature_importances_by_gain,
              feature_importances_by_weight,
              feature_importances_by_cover,
          ]
         ]
labels = ["permutation", "shap", "moris", "gain", "weight", "cover"]
colors = ["red", "blue", "green", "yellow", "pink", "grey"]
In [81]:
show_feature_saturation(scores, labels=labels, colors=colors)
/home/a_kulesh/anaconda3/envs/myenv/lib/python3.7/site-packages/plotly/tools.py:465: DeprecationWarning:

plotly.tools.make_subplots is deprecated, please use plotly.subplots.make_subplots instead


Sources